home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / DCLAP 6d / dclap6d / SeqPups / appsrc / autoseq.src / CTrace.H < prev    next >
Text File  |  1996-07-05  |  12KB  |  427 lines
  1. //    ============================================================================
  2. //    CTrace.H                                                        80 columns
  3. //    Reece Hart    (reece@ibc.wustl.edu)                                tab=4 spaces
  4. //    Washington University School of Medicine, St. Louis, Missouri
  5. //    This source is hereby released to the public domain.  Bug reports, code
  6. //    contributions, and suggestions are appreciated (to email address above).
  7. //    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -
  8. //    CTrace is a C++ template designed to store a sequence of numbers and perform
  9. //    simple operations on it.  It was originally designed for chromatogram data
  10. //    from automated DNA sequencing projects.
  11. //
  12. //    It was designed as a template for several reasons, but the most practical
  13. //    advantage is that we can generate a CTrace of doubles for the derivatives,
  14. //    even in cases where the original traces are shorts (for example).  One may
  15. //    then take the derivative of that (ie. the second derivative) by an identical
  16. //    call.  The behavior of using CTrace to store anything other than numbers is
  17. //    undefined at best, and at worst may subject you to immense ridicule from
  18. //    your peers (unless it works to solve the solution to the universe.)
  19. //
  20. //    CTrace maintains a pointer to an array of the parameterized type.  There are
  21. //    three ways to get data into this class:
  22. //        1. read the data yourself and pass the pointer to CTrace; don't forget
  23. //            to use the Size(ulong) method to specify the size. (This is space
  24. //            that you have allocated; do not pass temporary space to CTrace.)
  25. //            You should then call CalculateStats to complete the installation.
  26. //        2. open a file which contains a block dump of the specified type and
  27. //            pass the FILE* and size to read to ReadTrace. (CTrace allocates the
  28. //            space.)  By block dump, I mean something written with the equivalent
  29. //            of fwrite(YourDataPtr,ElementSize,#ofElements,FILE*).
  30. //        3. open an ifstream and pass it and the size to ReadAsText. (CTrace
  31. //            will allocate the space.)
  32. //    Data can also be written using the WriteTrace and WriteAsText methods.
  33. //    
  34. //    NOTES
  35. //    !    The destructor /always/ attempts to free space pointed to by the trace
  36. //        data member.  Set the trace pointer to NULL with Trace(NULL) if you
  37. //        wish to assume responsibility for disposal of the space.
  38. //    +    Statistics are automatically recalculated when necessary (ie. after
  39. //        setting the left clipping).
  40. //    *    Derivatives of an entire trace contain n-2 points; however, derivatives
  41. //        of clipped traces contain the full trace size.
  42. //
  43. //    SLATED IMPROVEMENTS
  44. //    *    better error mechanism (current method dirties namespace).
  45. //            - exceptions?
  46. //            - conserve global namespace by hiding error codes
  47. //    *    should ensure that peak list is empty before picking
  48. //    *    stats dirty flag in lieu of automatic recalc? (maybe faster...)
  49. //    *    derivative size inconsistency should be fixed.  If uncalculatable deriv
  50. //        values are filled, then it's wrong; but, we're forced to unnecissarily
  51. //        toss deriv values in the clipped trace just to maintain consistency.
  52. //        peakpicking may have to be modified if deriv numbering is changed.
  53. //
  54. //    MODIFICATION HISTORY
  55. //    93.06.30    Reece    Original coding
  56. //    93.07.29-31    Reece    Added Derivative, ReadAsText, WriteAsText methods
  57. //                        Converted to template, *AsText methods now use fstreams,
  58. //                        [] operator for trace access (returns lvalue).
  59. //    93.11.11    Reece    First release
  60. //    ========================================|===================================
  61.  
  62. #ifndef _H_CTrace                            // include this file only once
  63. #define _H_CTrace
  64.  
  65. #include    <stdio.h>
  66. #include    <iostream.h>
  67. #include    <fstream.h>
  68. #include    <math.h>
  69. #include    "CSequence.H"
  70. #include    "CPeakList.H"
  71. #include    "DNA.H"
  72. #include    "RInclude.H"                    // homegrown definitions
  73. #include    "RInlines.H"
  74. #include    "Definitions.H"
  75.  
  76. enum CTError
  77.     {
  78.     noError,                                // no error occurred
  79.     traceEmpty,                                // no trace to do that on
  80.     traceNotEmpty,                            // can't overwrite trace
  81.     badFile,                                // bogus file
  82.     memError,                                // memory error occurred
  83.     ioError,                                // i/o error occurred
  84.     rangeError,                                // something's out of bounds
  85.     analysisError                            // couldn't analyze peaks
  86.     };
  87.  
  88.  
  89. template<class T>
  90. class    CTrace
  91.     {
  92.     //private:                                // see note with Version(), below
  93.     //static vrsn    version;
  94.  
  95.     //
  96.     // Instance variables
  97.     //
  98.     CTError        error_flag;                    // error flag
  99.     T*            trace;                        // pointer to trace
  100.     size_t        size;                        // number of trace points
  101.     double        scale;                        // scale of trace
  102.     double        mean;                        // mean trace value
  103.     T            max;                        // max trace value
  104.     T            min;                        // min trace value
  105.     double        variance;                    // variance of trace values
  106.     T            baseline;                    // baseline correction
  107.     CPeakList    peaks;                        // list of peaks
  108.     tracepos    leftCutoff;                    // ignore indices before leftCutoff
  109.     tracepos    rightCutoff;                //   and after rightCutoff (0 based)
  110.     
  111.     //
  112.     // Methods
  113.     //
  114.     public:
  115.     //vrsn&        Version();                    // return class version (broken)
  116.                                             // there's a prize for figuring out
  117.                                             // why this works for other classes
  118.                                             // but not here...
  119.     CTError        Error(void);                // get and ...
  120. //    inline void    Error(CTError new_error=noError); //   set/clear error
  121.     inline void    Error(CTError new_error); //   set/clear error
  122.  
  123. //                CTrace(size_t sz=0);        // constructor, allocate space
  124.                 CTrace(size_t sz);            // constructor, allocate space
  125.                 ~CTrace(void);                // destructor
  126.  
  127.     CTError        AllocateTrace(size_t);        // allocate trace, update size, etc
  128.     inline void    DeallocateTrace(void);        // deallocate & update size, etc.
  129.  
  130.     inline size_t    Size(void);                // get and ...
  131.     inline void        Size(size_t new_size);    //   set trace size
  132.     inline T*        Trace(void);            // get and ...
  133.     inline void        Trace(T* tp);            //   set trace pointer (see NOTES)
  134.     inline CPeakList& Peaks(void);            // return reference to peak list
  135.  
  136.     inline T        Min(void);                // access functions for min
  137.     inline T        Max(void);                //   max,
  138.     inline double    Mean(void);                //   mean, and
  139.     inline double    Variance(void);            //   variance of trace values
  140.     CTError            CalculateStats(void);    // compute min,max,mean,variance
  141.  
  142.     inline tracepos    LeftCutoff(void);        // get and ...
  143.     inline void        LeftCutoff(tracepos t);    //   set left cutoff
  144.  
  145.     inline tracepos    RightCutoff(void);        // get and ...
  146.     inline void        RightCutoff(tracepos t);//   set right cutoff
  147.  
  148.     inline T&    operator[](tracepos index);    // trace value by index
  149.  
  150.     inline void        Baseline(T bl);            // get and ...
  151.     inline T        Baseline(void);            //   set baseline
  152.     inline double    Scale(void);            // get trace scale
  153.     void            Scale(double scale);    // scale trace; ie. 0.5 = scale 50%
  154.     void            Translate(T bl);        // add bl to each element
  155.  
  156.     CTError        Smooth(void);
  157.                 // Smooths by using a weighted average of the 3 points about a
  158.                 // particular index, with the special case of the endpoints
  159.                 // handled by throwing out the third point and normalizing the
  160.                 // weighting coefficients.
  161.  
  162.     CTError        Derivative(CTrace<double>** deriv);
  163.                 // Compute the derivative of the trace and return it in a
  164.                 // CTrace<double>.  deriv will be allocated automatically.
  165.                 // If deriv can be allocated, the derivative will be generated
  166.                 // and noError returned; otherwise, an appropriate error is
  167.                 // returned.
  168.  
  169.     CTError        PickPeakIndices(
  170.                     T        MinPeakHeight,
  171.                     double    ZeroThreshold,
  172.                     CSequence<tracepos>** peaks);
  173.                 // Picks local maxima>MinPeakHeight by a concave down method.
  174.                 // Returns in a CSequence<tracepos> all indices, i, for which
  175.                 // 1) the derivative at point i is within +/- ZeroThreshold or
  176.                 // 2) derivative at the point i-1 is positive and the derivative
  177.                 // at i+1 is negative (note that this implies the second
  178.                 // derivative is negative).
  179.  
  180.     CTError        PickPeaks(
  181.                     T        MinPeakHeight,
  182.                     double    ZeroThreshold);
  183.                 // Calls PickPeakIndices to get a sequence of tracepos.  This
  184.                 // sequence is converted into a sequence of peak records and
  185.                 // the bounds and width of the peak are computed.  The arguments
  186.                 // have the same meaning as those for PickPeaks.
  187.  
  188.     CTError        PeakBounds(
  189.                     tracepos    CenterOfSearch,
  190.                     T            Elevation,
  191.                     double&        LeftIntersection,
  192.                     double&        RightIntersection,
  193. //                    tracepos    MaxWidth = 0);
  194.                     tracepos    MaxWidth);
  195.                 // Determines the left and right bounds of a concave-down peak
  196.                 // at Elevation by searching laterally no more than MaxWidth on
  197.                 // either side of CenterOfSearch (ie. [COS-MW,COS+MW]).  If no
  198.                 // data point exactly equals Elevation, the approximate bound is
  199.                 // determined by linear interpolation.  The resulting left and
  200.                 // right bounds are returned in Left~ & RightIntersection.
  201.                 // MaxWidth=0 (default) specifies no limits (ie. search is
  202.                 // exhaustive, limited only by limits of trace).  If the search
  203.                 // bounds [COS-MW,COS+MW] exceed the bounds of the trace, they
  204.                 // are automatically truncated to the nearest limit.
  205.  
  206.  
  207.     // I/O
  208.     CTError        WriteAsText(ofstream& os);
  209.                 // Writes trace from a stream as a return-delimited list of
  210.                 // points.  All numerical data types are supported.
  211.     CTError        ReadAsText(ifstream& is, size_t);
  212.                 // Read analog of ReadAsText.
  213.  
  214.     CTError        WriteTrace(FILE*);
  215.                 // Writes trace as a block of numbers.  The good news is that
  216.                 // it's extremely fast; the bad news is that the file is
  217.                 // essentially a dump of the internal representation of the
  218.                 // numerical array and is difficult (for humans) to interpret.
  219.     CTError        ReadTrace(FILE*, size_t);
  220.                 // The read analog to the above.
  221.  
  222.     friend ostream& operator<<(ostream& os, CTrace<T>& t);
  223.                 // Writes a user-friendly summary of CTrace data members and
  224.                 // the list of picked peaks if there are any.
  225.     };
  226.  
  227. //
  228. //    INLINE FUNCTION DEFINITIONS
  229. //    ===========================
  230. //
  231.  
  232. //template<class T>
  233. //inline
  234. //vrsn
  235. //CTrace<T>::Version()
  236. //    {
  237. //    return version;
  238. //    }
  239.  
  240. template<class T>
  241. inline
  242. void
  243. CTrace<T>::DeallocateTrace(void)
  244.     {
  245.     delete trace; size = 0; min = max = 0; mean = 0;
  246.     }
  247.  
  248. template<class T>
  249. inline
  250. CTError
  251. CTrace<T>::Error(void)
  252.     {
  253.     return error_flag;
  254.     }
  255.  
  256. template<class T>
  257. inline
  258. void
  259. CTrace<T>::Error(CTError new_error)
  260.     {
  261.     error_flag = new_error;
  262.     }
  263.  
  264. template<class T>
  265. inline
  266. double
  267. CTrace<T>::Mean(void)
  268.     {
  269.     return mean;
  270.     }
  271.  
  272. template<class T>
  273. inline
  274. double
  275. CTrace<T>::Variance(void)
  276.     {
  277.     return variance;
  278.     }
  279.  
  280. template<class T>
  281. inline
  282. T
  283. CTrace<T>::    Max(void)
  284.     {
  285.     return max;
  286.     }
  287.  
  288. template<class T>
  289. inline
  290. T
  291. CTrace<T>::    Min(void)
  292.     {
  293.     return min;
  294.     }
  295.  
  296. template<class T>
  297. inline
  298. T&
  299. CTrace<T>::operator[](tracepos index)
  300.     {
  301.     return trace[index];
  302.     }
  303.  
  304. template<class T>
  305. inline
  306. void
  307. CTrace<T>::LeftCutoff(tracepos t)
  308.     {
  309.     if ((t>=0) && (t<=size-1))
  310.         {
  311.         leftCutoff = t;
  312.         CalculateStats();
  313.         }
  314.     }
  315.  
  316. template<class T>
  317. inline
  318. void
  319. CTrace<T>::RightCutoff(tracepos t)
  320.     {
  321.     if ((t>=0) && (t<=size-1))
  322.         {
  323.         rightCutoff = t;
  324.         CalculateStats();
  325.         }
  326.     }
  327.  
  328. template<class T>
  329. inline
  330. CPeakList&
  331. CTrace<T>::Peaks(void)
  332.     {
  333.     return peaks;
  334.     }
  335.  
  336. template<class T>
  337. inline
  338. size_t
  339. CTrace<T>::Size(void)
  340.     {
  341.     return size;
  342.     }
  343.  
  344. template<class T>
  345. inline
  346. void
  347. CTrace<T>::Size(size_t new_size)
  348.     {
  349.     size = new_size;
  350.     }
  351.  
  352. template<class T>
  353. inline
  354. T*
  355. CTrace<T>::Trace(void)
  356.     {
  357.     return trace;
  358.     }
  359.  
  360. template<class T>
  361. inline
  362. void
  363. CTrace<T>::Trace(T* tp)
  364.     {
  365.     trace = tp;
  366.     }
  367.  
  368. template<class T>
  369. inline
  370. double
  371. CTrace<T>::Scale(void)
  372.     {
  373.     return scale;
  374.     }
  375.  
  376. template<class T>
  377. inline
  378. void
  379. CTrace<T>::Baseline(T bl)
  380.     {
  381.     baseline = bl;
  382.     }
  383.  
  384. template<class T>
  385. inline
  386. T
  387. CTrace<T>::    Baseline(void)
  388.     {
  389.     return baseline;
  390.     }
  391.  
  392. template<class T>
  393. inline
  394. tracepos
  395. CTrace<T>::LeftCutoff(void)
  396.     {
  397.     return leftCutoff;
  398.     }
  399.  
  400. template<class T>
  401. inline
  402. tracepos
  403. CTrace<T>::RightCutoff(void)
  404.     {
  405.     return rightCutoff;
  406.     }
  407.  
  408. template<class T>
  409. ostream&
  410. operator<<(ostream& os, CTrace<T>& t)
  411.     {
  412.     tracepos length = t.rightCutoff-t.leftCutoff+1;
  413.     os    << "trace size:\t"
  414.             << t.Size() << endl
  415.         << "Window:\t\t"
  416.             << "[" << t.leftCutoff << "," << t.rightCutoff << "]"
  417.             << "  (" << length << " point" << Plural(length) << ")" << endl
  418.         << "min/max/mean/var: "
  419.             << t.Min() << "/"
  420.             << t.Max() << "/"
  421.             << t.Mean() << "/"
  422.             << t.Variance() << endl;
  423.     return os;
  424.     }
  425.  
  426. #endif                                        // conditional inclusion
  427.